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Abstract. A direct three dimensional FIT reconstruction algorithm based on 
complex geometrical optics solutions and a nonlinear scattering transform is presented 
and implemented for spherically symmetric conductivity distributions. The scattering 
transform is computed both with a Born approximation and from the forward 
problem for purposes of comparison. Reconstructions are computed for several test 
problems. A connection to Calderon's linear reconstruction algorithm is established, 
and reconstructions using both methods are compared. 



Direct numerical reconstruction of conductivities in 3-D 



2 



1. Introduction 

The reconstruction of conductivity distributions in two or three dimensions from 
measurements of the current density-to- voltage map is known as electrical impedance 
tomography, or EIT, and has applications in medical imaging, nondestructive testing, 
and geophysics. For the 3-D bounded domain considered here, medical applications 
include head imaging and the detection of breast tumors. See, for example, [Hol05] for 
a survey of clinical applications of EIT. In this work, we consider a bounded domain 
in and present a direct reconstruction algorithm and its numerical implementation 
on the unit sphere. The theoretical foundation of the method dates back more than 20 
years to a series of papers by Sylvester- Uhlmann [SU87], Novikov [Nov88], Nachman- 
Sylvester-Uhlmann [NSU88] and Nachman [Nac88]. The algorithm makes use of complex 
geometrical optics (CGO) solutions to the Schrodinger equation and uses the inverse 
scattering method. This is described in detail in section 2 of this paper. 

The inverse conductivity problem was first formulated mathematically by Calderon 
[CalSO] as follows. Let Q C R",n > 3 be a simply connected, bounded domain with 
smooth boundary dfl, and let 7 G L°°{fl) denote the conductivity distribution. Assume 
there exists C > such that for x G Q, < 7(2;) < C . The electric potential u 
arising from the application of a known voltage to the boundary of Vt is modeled by the 
generalized Laplace equation with Dirichlet boundary condition 



Thus A-j, represents static electrical boundary measurements: it maps an applied voltage 
distribution on the boundary to the resulting current flux through the boundary. 
Calderon [CalSO] posed the question of whether the conductivity 7 is uniquely 
determined by the Dirichlet-to-Neumann map, and if so, how to reconstruct the 
conductivity. He gave an affirmative answer to the uniqueness question for the linearized 
problem and gave a reconstruction algorithm for that case. His algorithm is described 
in section 2.3 of this paper. 

The uniqueness question for 7 G L°°{Vt) is still open in R'^, but has been solved 
recently by Astala and Paivarinta [AP06b] for a bounded domain in R^, sharpening 
the previous results due to Nachman [Nac96] in which 7 G p > 1, and Brown 

and Uhlmann [BU97] in which 7 G p > 2. In three dimensions the uniqueness 

problem was solved for smooth conductivities in [SU87] . At the time of this publication, 
in R^ the uniqueness results with lowest regularity are [BT03] with 7 G 
p>2n and [PPU03] with 7 G W^/^'^{n). 

Most existing 3-D EIT reconstruction algorithms are linear or iterative, minimizing 
a functional that describes the nearness of the predicted voltages to the measured 
data in a given norm with one or more regularization terms. In contrast, the 
algorithm presented here is direct and fully nonlinear. It is similar to the 2-D D-bar 



V ■ 7VM = in f2, M = / on 9^2. 
The Dirichlet-to-Neumann map A^ is defined by 



(1) 




(2) 
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algorthims based on the works [Nac96] and [BU97], which were first implemented in 
[SMIOO, SMIOl, Knu03, MS03]. In these initial works the Born approximation to the 
CGO solutions is used in the computation of the scattering transform. It was used 
successfully on experimental tank data in, for example, [IMNS04, EM09] and human 
chest data in [IMNS06, DM10]. This inspired the approach in section 2.2 of this paper in 
which the Born approximation is used in the 3-D direct algorithm. For further reading 
on 2-D D-bar algorithms, the reader is referred to [KLMS07] in which the application 
to discontinuous conductivity distributions is specifically addressed, and [KLMS09] 
in which a rigorous regularization framework is established using the full scattering 
transform. Calderon's method has also recently been used for the reconstruction from 
experimental data in both 2-D [BM08] and 3-D [BTJIS08]. 

In this work, we assume the conductivity 7 G C^{Q), we take Q to be the unit sphere 
in R^, and assume 7 = 1 near dQ. The smoothness assumption on 7 is necessary, but 
the other assumptions are made mainly for simplicity in the numerical computations. 
We stress in particular that the theory is valid in more complex geometries. The study 
of the effects of noise in the data is not in the scope of this paper, but rather is left for 
future work. 

The outhne of the paper is as follows. In section 2.2 we describe a direct 
reconstruction algorithm with a linearizing assumption tantamount to a Born 
approximation. That approach is referred to as the t'^'P approach, consistent with the 
notation used in the 2-D D-bar algorithms. An explicit connection to the linearized 
method of Calderon is established in section 2.3. The reconstruction of the conductivity 
in the 2-D D-bar method described in the works above is achieved by taking a small 
frequency limit in a D-bar equation for the CGO solutions to directly obtain 7(2;). In 
contrast, here we have to take a high complex frequency limit. A D-bar equation for 
the 3-D problem is utilized in [CKS06], resulting in a promising, but more complicated 
approach than the one studied here. The numerical implementation of that approach 
is left for future work. In section 3 we consider the case of spherically symmetric 
conductivities and show symmetry properties in the scattering transform. We also show 
how the Dirichlet-to-Neumann map can be represented and approximated in that case. 
Details on the numerical implementation are found in section 4. Numerical examples 
are found in section 5. 

2. The reconstruction methods 

In this section we describe the theoretical reconstruction method, the t"'''' approach, and 
the relationship to Calderon's linearized method. 

2.1. The nonlinear reconstruction method 

The method was developed in [SU87, NovSS, NSU88, NacSS]; here we provide a brief 
outline. The reader is referred to [Nac88] for rigorous proofs. The equations closely 
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parallel those of the 2-D problem (note that [Nac88] precedes that work), and so readers 
familiar with that case will recognize the notation and functions involved. 

The initial step is to transform the conductivity equation into a Schrodinger 
equation. Indeed, if u satisfies (1) then v = satisfies 

(-A + = in Q with q^^—. (3) 
Note that q = near dfl. The Dirichlet-to-Neumann map for equation (3) is defined by 

/. , dv 

where now v satisfies (3) with v\qq — f . In general the maps and are related by 

The assumption that 7 = 1 in a neighborhood of dfl simplifies (4) to Ag = A^. 

To define the CGO solutions, introduce a complex frequency parameter C £ and 
define the set 

V = {C3\{0}:C-C = 0}. (5) 

Then e*^'^ is harmonic in if and only if C e V. For ^ e R^, introduce the subset of V 
given by 

H = {C e V: (e + 0' = 0}. (6) 

Note that ( ■ ( — (C + C)^ — gives an explicit characterization of in terms of 
an auxihary vector ^-^ e R^ with • ^ = 0. Indeed suppose (r, G R^. Then 
( — (u + iQ E if and only if 

Cr =-e/2 + e^ (7) 

Since g = in a neighborhood of dQ,, one can extend q — into R^ \ Q. The CGO 
solutions ip{x, to the Schrodinger equation solve 

(-A + g(x))^(x,C) = 0, xGR^ CeV, (8) 

and behave like e*^ '' for \Q large. More precisely, define 

Then fi — 1 approaches zero in a certain sense as cither \x\ or \(\ tends to infinity, see 
[SU87, Nac88]. Note that ■ip{x,C) grows exponentially for x ■ lm( < 0. The function /j, 
satisfies 

(-A - 2< ■ V + C) = in R3. (9) 

Denote by the Faddeev Green's function defined by 

G,ix)^e'-^g,ix), where 5c(^) = ^Y^^^C, (10) 
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where the integral is understood in the sense of the Fourier transform defined on the 
space of tempered distributions. The functions G^, g^^ are fundamental solutions of the 
Laplace equation and conjugate Laplace equation respectively, i.e. 

A^gr^ = -So and AG^ = -So- (11) 

Then (9) is equivalent to the Faddeev-Lippmann-Schwinger equation 

(I + g<:*(q-))l^^linR\ (12) 

Estimates for the operator g^^* for large ( ([SU87]) and small C ([CKS06]) give the 
existence and uniqueness of /i (and therefore ip) for any sufficiently large or small C ^ ^■ 
The key intermediate object in the reconstruction method is the so-called non- 
physical scattering transform of the potential q defined for ^ e and sufficiently large 
or small ^ G V by 

t(e,C)= / e-'^<^+^^^l;(x,Oq(x)dx. (13) 
Jn 

Integrating by parts and assuming that ( eV^ we find that 

t(e, C) - / e— «+^)(A, - AoMx, C)da{x). (14) 

Thus we require ■0|an iii order to compute the scattering transform from the Dirichlet- 
to-Neumann map. It turns out that 'ipldci satisfies a uniquely solvable Fredholm integral 
equation of the second kind on dfl [Nov88, NacSS], namely, 

7lj{x,C)+ [ G^{x-x){Ag-Ao)7/j{x,C)da{x)^e'''<, x e dQ. (15) 
Jan 

Note from (13) and the fact that ijj ~ e'"''^ that from the scattering transform one 
can compute the Fourier transform q of the potential by taking the large frequency limit 

hm t{^,0 = m- (16) 

|C|->oo 

Summary of the reconstruction method: 

(i) Solve the boundary integral equation (15) for ■0|an- 

(ii) Compute t(e, C) for ^ e R^, C e H by (14). 

(iii) Compute q{^) from (16) . 

(iv) Compute q by inverting the Fourier transform. 

(v) Compute 7 by solving -A^ + q^^ = in Q, ^y^\an = 1- 

We stress that the ill-posedness of the inverse problem is in this algorithm isolated 
in the first step. 

2.2. The reconstruction method using t""^ 

Inspired by the t"'''' approximation in the 2-D D-bar method, an analogous approach 
can be taken in 3-D. Approximating il'{x, Q on the boundary by its asymptotic behavior 
gia; C ehminates the need for the ill-posed first step. We define for ^ e R^, C £ 

t^^'(^,C)=/ e-'^<^+^\Ag-Ao)e'^<da{x). (17) 
Jan 
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This approximation is tantamount to a linearization of the first step in the reconstruction 
algorithm above around 7 = 1. Using for t in (16) gives the following simple 
reconstruction algorithm: 

(i) Compute t--(e, C) for ^ e R^ C e H by (17). 

(ii) Compute 

r^(0= lim t^-(e,C) (18) 

ICKoo 

and then q"'''^ by inverting the Fourier transform. 

(iii) Compute 7'='''' by solving —A^/y^ + g^^p^T^ = in 1] y^T^\dn = 1- 

It is not guarenteed from the theory that the limit in (18) is well-defined. In our 
numerical simulations we will compute t'"'^((^,^) for a fixed but large value of (. This 
will numerically define q'"''^{^)- 

2.3. Calderon's linearized reconstruction method 

Several properties of t™"" can be established from an analysis comparing this approach 
to that of Calderon. In [KLMS07] a connection was established between the 2-D D- 
bar method based on the global uniqueness proof by Nachman [Nac96] and Calderon's 
linearized reconstruction method. 

Define a function u'''^^{x, Q as the unique solution to the boundary value problem 

V • 7VM'"'^(a;, C) = 0, x e Q, C e 

Integration by parts in equation (17) results in a formula for t"'''" defined in terms of 7 
in the interior 

t-p(^, () = /" (7 - l)VM^''^(a;, C) • Ve-'^-(«+'^)(ix. (19) 
Jo. 

Write li'^'^p = e^^< + 5u for 5u e Hl{Vt). Then 5u satisfies 

V • {^V5u) = - V • ((7 - l)Ve^^-^), (20) 
and one can estimate 

¥u\\H^m < C\\{^ - l)Ve''-%Hn) < ICIII7 - l||Loo(n)el«l^ (21) 
where R is the radius of the smallest ball containing fl. From (19) we then get 

t-p(^, C) = / (7 - l)V(e^^-^ + 5u) ■ Ve-'^-^^+^^dx 
Jn 

= / (7 - 1) Ve^^< • Ve-^^-(«+^)cia; + R{^, Q 
Jn 

= (e-C) / (7-l)e-'^'^(ix + i?(e,C), 
Jn 

where the remainder term 

i?(^, = /" (7 - i)VSu ■ Ve-'^'-^^+^^dx. 
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Since + — — we have — — 2^ ■ ( and hence 

t-'i^X) = + ^(e,C). (22) 

The remainder is 0(|C|) for ( small, which can be seen from (21). This fact suggests 
that we use the minimal C e V^, that is 

a = withO-e = 0, |C/| = y- 

Moreover, with this particular choice we can divide in (22) by |^|^ as the following 
proposition shows. 

Proposition 2.1. Suppose 7 e L°°(Q). Then 

it-ua)i = (^(i^r) 

for small \^\. 

Proof. Note that = |CP/2. Since maps constant functions to zero and has its 
range inside the space of mean free functions in H~^/'^{dQ), we have that for small |^| 

I (e-"-(5+C«) _ i)(A^ _ A^)(e«<e - l)da{x) 
Jan 

and hence the particular form of Q gives the conclusion. □ 



|t^'''(e,a)i 



With the particular choice ( = Q in (22) we now neglect the term R{^, and 
divide by — which gives 

Introduce Xb{0^ characteristic function on the ball |^| < B. With this function 
we remove high frequencies and invert the Fourier transform. This results in a linear 
reconstruction algorithm 

This formula is equivalent to the second inversion formula obtained by Calderon [Cal80, 
p. 72]. 

In summary the linear reconstruction algorithm consists of two steps: 

(i) Compute t-''(e,a) by (17). 

(ii) Compute the reconstruction by (23). 

This method is truly a linearization of the nonlinear reconstruction method outlined 
in section 2.1. As explained above t™^ is a linearization of the first step on page 5. 
Moreover, the computation of the quantity 



(27r)= 

linearizes the step q ^ Finally, linearizing the square function gives (23). 
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3. The case of a spherically symmetric conductivity 

As a test problem it is of special interest to consider spherically symmetric conductivities 
in the unit sphere. In this case the scattering transform has certain symmetry properties 
described below. Moreover, the Dirichlet-to-Neumann map is described explicitly in 
terms of eigenvalues and eigenf unctions, which in this case are the spherical harmonics. 
These properties will be derived in this section. 



3.1. Symmetry in the scattering transform 

The Fourier transform of a spherically symmetric function is spherically symmetric itself. 
For the scattering transforms t and we have similar prperties. In the following we 
will tacitly assume that ( is either small or large such that t(^, () is well-defined. 

Proposition 3.1. Let R e SO{2) be arbitrary, and suppose q{x) — q{Rx) for x E Q. 
Then 

t(e, C) = t(i?e, RC), t-(e, C) = t-^(i?e, RC) (24) 

In particular, 

m, Ci) = t(e, C2), t-(e, Ci) = t-^{^, C2) (25) 

for all (1X2 e V^. 

Proof. We will prove the result for t only; for f""^ the reasoning is similar. Let 
R e SO{2). By the uniqueness of the CGO solutions, the rotational invariance of the 
Laplace operator, and the symmetry in q we have ip{xX) — i^iR^, RQ- Consider the 
integral (13) and make the change of variables R^y — x : 



t(e,C)= / e-'^<^+^^ij{xX)q{x)dx 

e-'''^y<^^MR^yX)<i{R^y)d{R^y) 

'RCl 

= f e-'y^^^+^my,RQq{y)dy 
Jn 

^t{R^,RC). 

To prove (25) fix ^ e and take (1,(2^ with \(^\ = {(^l- Then 

= Re(0) + ilm(O), Re(0) = + • C = 0, 

lm(0) • e = lm(0) • = 0. 

Define a linear transformation R by R^ = ^, R^^ = -^(C x Ci^) = C x ^2"- a 
consequence R e SO{2) and (25) follows from (24). □ 



The Fourier transform of a real and even function is real itself. For the scattering 
transform we have the following equivalent property: 
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Proposition 3.2. Suppose q{x) is real and even. Then for ^ e R^, C ^ H 



t(e,C) = t(e,C), t-(e,C) = t-(e,C)- (26) 

Proof. We will again only show the properties for t. From the uniqueness of the CGO 
solutions it follows that if q is even then ij,{—x,C) = A*(a;, — C)- Moreover, if q is real, 
then ii{x, Q = fi{x, —(). Hence, if q is both even and real then 

WZ)- I e-«g(x)//(a;,C)o?x= / e-^^-^y^q{x)yi{x,-Qdx 

= f e-'y<q{y)ii{yX)dy 
Jn 

= t(e,C)- 



□ 



We now have a corollary for spherically symmetric potentials: 
Corollary 3.3. Suppose q is spherically symmetric. Then 



t(e,C) = t(e,C), eeR', CeVe 

Proof There exists R e SO (2) such that 

i?(0=e, R{0 = C 

and hence from (24) we have t(,^, () — t(^, Q. Equation (26) now implies the result for 
t. For the result follows similarly. □ 

3.2. Eigenfunctions and eigenvalues for the Dirichlet-to- Neumann map 

We use the same ideas for the computation of eigenvalues for the 3-D problem that were 
used for the 2-D problem in [SMIOO]. 

Proposition 3.4. Let D be the unit disk and suppose ^{x) is spherically symmetric. 
Then the eigenfunctions of are the spherical harmonics Yf^. 

Proof. When 7 is spherically symmetric it follows from separation of variables that the 
solution to V • "y'Vuim — with uim\dD — Yi"^ 

uim^Ri{r)Yr{e,<j)), (27) 

where Ri{r) solves an Euler type equation. Thus 

dRi 



dr 



Yr{eA)^kY{^{OA). (28) 

r=l 

□ 



Note that A is independent of m since Ri is independent of m. 
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3.3. Approximation of Eigenvalues and Eigenfunctions of the Dirichlet-to-Neumann 
Map 

Next we will consider how to approximate the eigenvalues for the special case of a 
constant conductivity 7 = 1. The particular form of Ri gives the following result. 

Proposition 3.5. The eigenvalues of Ai are given by Xi = I. 

In the case of a piecewisc constant radially symmetric conductivity the eigenvalues 
can be computed recursively. Suppose = ro < ri < r2 < . . .rjv_i < r^ — 1 and for 
j^l,2,...,N 

l{x)=-fj>0, \x\e[rj_^,rj]. (29) 

Proposition 3.6. Suppose 7 is given by (29). Then the eigenvalues of are given by 

2/ + 1 

Ao = 0, Xi^l- — — ^ , / > 



where = ^ ^jnm£i±li ^ith pi = 1, p. = , A = ¥ and w. = rT 



{21+1) 



- IAn -{1 + 1)Bn =l-{2l + 1)Bn (30) 

dD 



Proof. Since is a constant, Aq = 0. The solution to V • ^Vuim — 0, ui^nldo — Y/^, 
is given by (27) with Ri(r) = Ajr'- + Bjr~^^+'^^ for r^-i < r < rj, j = 1,...,A^. 
The coefficients Aj and Bj are determined by matching the Dirichlet and Neumann 
conditions at the r^, j = 1, . . . , — 1. The outermost Dirichlet condition (at r — 1) 
gives 1 — A]si -\- B]si which leads to the following eigenvalue expression: 

Xi = l 

or 

Moreover, by induction it follows that Aj = BjCj^i for j = 2, . . . , N. Again using 
the Dirichlet condition from the boundary, 1 = Aj^ + B^ wc get B]\f = (Cjv-i + 1)""^ 
which leads to the expression of the eigenvalue as stated in the theorem. □ 

By [SCII91] if conductivities jl and 7?/ are such that jiir) < 7[/(r) for all r, then 
the eigenvalues Xf and X^ of their corresponding Dirichlct-to-Neumann maps satisfy 
Af < Xf . This gives a means for finding lower and upper bounds on the eigenvalues of 
a smooth function by finding the eigenvalues of piecewise constant function, 7^ and ju 
that satisfy 7L(r) < j{r) < ■ju{'>^)- 



4. Implementation details 

4.1. Numerical method for computing the scattering transform t(^,C) 

We compute the scattering transform t(^,C) from the definition (13) as a comparison 
to the f'^p approximation and to study the reconstructions from an accurate scattering 
transform. The computation requires that we solve the Lippmann-Schwinger equation 
(12) for ii{x, Q. Hence we require 

• A method of computation for the Faddeev Green's function in three dimensions 
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• A numerical method for the solution of (12) 

• Numerical quadrature for computing t(^, Q from (13) 



We describe each of these in turn. 

4.1.1. Computation of the Faddeev Green's function The Faddeev Green's function was 
defined in equations (10). The effect of scaling and rotation of ( on G(; was analyzed in 
[CKS06], and it was shown that when ( satisfying ^ • ^ = is decomposed in the form 

(^K{k± + ik), (31) 
where k±, k e R^, — \k\ — l,k ■ k± — 0, and |C| = V^n, then 

g<;{x) = K'^'^gk^+ikiKx). (32) 
Furthermore, if i? e SO{2) then 

gd^) = gnd^x). (33) 

Combining (32) and (33) yields the formula 

gdx) = Kge^+iedf^-Rx): (34) 

where R e SO (2) and the first and second column of R is k,k± respectively. This 
formula shows that it is sufficient to compute (7ei+ie2- 

To compute gei+ie2 we will use formula (6.4) of [New89] 

ger+ie2{x) = — / — ^ Ji{rVl - u'^)du, (35) 



4nr 4:71 J, ^\-u^ 

where J\ denotes the Bessel function of the first find of order one. Here r = and 
s = X ■ e-i = xl\x\ ■ 62. Since the function Ji{t)/t is continuous on the interval [0, oo) 
(in particular at t = 0), we will approximate the integral in (35) by a simple midpoint 
Riemann sum 

/■I p—ru+X2—ixi ^ p—ru(j)+X2—ix\ 

where N is the number of discretization points, h = {1 — s)/N and u{j) = s + {j — 
l/2)h, J = 1,2,. ..,7V. 

4- 1.2. The computation of complex geometrical optics Having computed the Faddeev 
Green's function we now turn to the numerical solution of the integral equation (12) for 
/x(-,C). We will use a method due to Vainikko [VaiOO] for solving Lippmann-Schwinger 
equations; see also [HohOl, KMS04] for implementations in different contexts. The main 
idea is to transform (12) to a multiperiodic integral equation in R^, which can be solved 
efficiently using FFT. 

Let Gp = {x I l^il < p}. Then by assumption supp(g) dVL d Gi. Extend the 
potential q and the Green's function gc^ to G2 such that 

^^^"10, xeG2\n, ^^^^^ " 1 0, x e 6-2 \ 1^, 
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and then extend and to as periodic functions in all variables with period equal 
to 4. Instead of (12) we consider the periodic integral equation 



f^''ix,C)+ glix-y)q^iy)fi^iyX)dy^l. (36) 

This equation is uniquely solvable since (12) is, and moreover one can show that on fl 
we have 

IJ,''{x,C) ^ li{x,C), xen. 
In order to solve (36) numerically define 

Z% = {jeZ^\ -N/2<jk<N/2, k = l,...,3} 
and the computational grid 

where h — A/N specifies the discretization fineness. Define the grid approximation 0jv 
of a continuous function e C(G2) by 

0iv(i/i) = (l)(jh) 

and the grid approximation g^ of g^ (which is smooth except for a singularity at the 
origin) for fixed ( by 



9N{jh) 



0, j = 

g'^ijh), otherwise. 

The convolution operator appearing in (36) 



K(f){x) = / gl{x - y)(t){y)dy 



is now discretized by trigonometric collocation, which, using the discrete Fourier 
transform J^at, gives 

KN(t)Nijh) = J'N^igN ■ 'Pn)- 

Here • denotes pointwise multiplication. In practice the discrete Fourier transform 
can be implemented efficiently using FFT (with proper zero-padding) in 0{N^ log{N)) 
arithmetic operations. The total discretization of (36) now reads 

This discrete linear system is solved numerically in matlab using the iterative algorithm 
GMRES [SS86], without setting up a matrix for the hnear map KN{qN-). 



4.1.3. The scattering transform Having computed the grid approximation /i^r it is 
straightforward to evalute t{^,() by using numerical integration in (13). In this 
implementation we have used a simple midpoint quadrature rule. 
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4-2. Numerical method for computing t""'' for spherically symmetric conductivities 

For the calculation of we expand e'^'^ in terms of spherical harmonics^ and e""'^'^^'*'^^ 
in terms of the spherical harmonics conjugates, 

oo / 

1=0 m=-l 
oo k 

k=0 n=—k 

Using these expansions leads to 

l,m k,n -^9^ 

In the special case of spherically symetric conductivities we can use the knowledge 
of the eigenvalues of the Dirichlet-to-Neumann maps, A-yYi"^[6, (j)) = XiYf^ to simplify 
the calculation of t''''''. In particular we get 



l,m,k,n -^^^ 

= E «^-(^'C)MC)(A.-A;) /" 

l,m,k,n -^^^ 

= E«^-(^'OMC)(A/-0 (38) 

l,m 

The last equality comes from the orthonormality of the spherical harmonics. Equation 
(38) can be easily calculated if the coefficients a^^ and bim are available. In this work, 
these coefficients were calculated with a software package called ''S2kiV which are C 
routines that can be accessed from Matlab. Detailed information can be found in 
[HRKM03]. 

4-3. Computation of the conductivity 

After taking the high frequency hmit in (16) and (18) we calculate the inverse Fourier 
transform to get q{x) and q^'^'^{x). The integral in the inverse Fourier transform is here 
computed numerically using a simple Riemann sum. To get the conductivity 7 we need 
to solve the boundary value problem A7^/^ = 57^/^ with 7^/^|an = 1- This was realized 
with the standard Green's function for the Laplace equation in three dimensions. Using 
symmetries reduces the problem to a single integral. 

4.4- Numerical implementation of Calderon's method 

Calderon's method based on (23) is simply implemented by evaluating the integral using 
numerical quadrature. 

I We use here the normahzed spherical harmonics given by Y{^{6,(f)) = Nf" PI" {cos 6) e^"^"^ where A^" 
are normahzation factors and are associated Legendre functions. 



Direct numerical reconstruction of conductivities in 3-D 



14 



5. Results 

5.1. Examples 

The conductivity distributions we will use in the examples are smooth, spherically 
symmetric and constant one near dQ. They are given by 

7(x) = {a^{\x\) + lf 




for — d < r < d 
otherwise 



(39) 



where < d < 1 is a parameter determining the support of ^. The parameter a regulates 
the amplitude of 7, which is largest at r = with amplitude (a+l)^. A similar function 
was used [SMIOO] as an example for the two dimensional problem. 

5.2. The scattering transform 

Let us fix (i = 0.9, a — 0.3. in (39). We are interested in the limit of t(^, C), C ^ R-^? C ^ 
V^, when \(^\ goes to infinity. For purpose of illustration we compute t{^,() fo^' fixed 
^ = (10,0,0) and varying ( E with 8 < |C| < 50. We use a discretization level 
in the algorithm corresponding to = 2®. In addition we compute t""''^ {C, , () by (37) 
using the first 30 eigenvalues of the Dirichlet-to-Neumann map. We truncate the sum 
of the spherical harmonics at / = 30, which means we use approximately the first 900 
spherical harmonics. As a benchmark we compute q{i). The results are shown in figure 
1. We know from Corollary 3.3 that t and are real and this is consistent with our 




Figure 1. t(^,C), t«"P(^,C) calculated for fixed ^ 
d = 0.9 and a = 0.3. 



(10,0,0) and varying \(^\. Here 



numerical results. The data verifies that for our example t(^,^) converges to q{^) as 
C ^ 00. We observe that is independent of the magnitude of C G Vj, until it diverges 



Direct numerical reconstruction of conductivities in 3-D 



15 



due to numerical instability. The same phenomena appears in other examples and with 
different values of ^. We believe that this phenomena has to do with the special class of 
spherically symmetric conductivities considered here. 

Next we compare t{^,() ^-iid t'"'p(^,(^) for different values of ^. For each ^ = 
s[l,0, 0], s e [0,50], we fix C e with \(\ = 50. We compute t(^,C) using a 
discretization level with N = 2^. t^""^ is computed with the parameters as above. As a 
benchmark we compute q{i)- The results are displayed in figure 2. The difference in q{^) 

1 



0.5 





-0.5 

-1 
-1 .5 

-2 

10 20 30 40 50 

Figure 2. Scattering data t, f'^P and q (qhat) with d = 0.9 and a = 0.3. For each ^, 
C € V{ is chosen such = 50. The Fourier transform q virtually coincides with t(^, (). 

and t(^, is very small. t'"'p({) is displayed only for < |^| < 32 since the calculation 
becomes numerically unstable and blows up for |^| > 32. One observes good agreement 
of all three curves for |,^| > 20. Close to |^| = the approximation t*""'" is close to zero 
and differs from the correct values. 

5.3. The reconstructions 

Evaluting the inverse Fourier transform of the numerically computed t'"'''(,^) and t(^, () 
gives two approximations of q{x) which are displayed in figure 3. The approximation 
calculated from t(^, ^) differs as expected only slightly from the actual value. The 
approximation q'"'^ of q calculated from t^^^ (and hence from the boundary data) is quite 
different from q. For x near the boundary the q'"'^{x) is quite accurate, but for x near 
zero there are large discrepancies, especially in the magnitude. Looking at the scattering 
data in figure 2, one sees two features most likely responsible for that difference. The 
first one is the differences in the values of t'"'p(^) for ^ close to zero compared to q{^). 
The second is the truncation of t'"'p(^) due to numerical instability for large ^ values. 
More details on the influence of the truncation of the scattering data are provided in 
section 5.4. 
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Figure 3. Left: Reconstructions of q{x) by taking the inverse Fourier transform of 
t(^, () and t™P(^, () for a = 0.3 and d = 0.9. Right: Reconstructions of 7 from t, -y'^'^p, 
and 7*PP compared to actual conductivity for a = 0.3 and d = 0.9. 7 from t nearly 
coincides with the actual conductivity. 



Also in figure 3 we display three reconstructions of the conductivity distribution. 
The first reconstruction of 7 is from t{^, Q. Since t{^, Q is computed from the forward 
problem, it may be expected that this reconstruction would be very close to the actual 
value, as it is. The second reconstruction is Y^^^i^) from and the third reconstruction 
^app jg fpQjj^ ^jjg linear method (23). Considering the relatively large difference in 
magnitude of q''''^{x), the reconstruction 7^^^^ is surprisingly good. Also 7^^^ is a fairly 
good reconstruction. A positive aspect in both reconstructions is that wc get 7 = 1 
close to the boundary. Moreover, the overall shape is also fairly well reconstructed. 

5.4- The influence of the truncation of the scattering data 

When we reconstructed q'"'^ and 7""'' we truncated the scattering data t""'' due to 
numerical instabilities. In this section we investigate the influence on the reconstructions 
of the truncation of the true scattering data t{C,,()- Figure 4 shows t(^,(^) and the 
reconstructions q{x) and 7(2;) for different truncations of t(,^,(^), namely at ,^ = i? 
for R = 15,25,50. We have chosen ( E with \(\ = 50. The actual potential and 
conductivity are almost identical to the curves corresponding to R — 50. It is evident 
that the amount of truncation of the scattering transform influences the reconstruction, 
and that a very poorly reconstructed q can still result in a good approximation of 7. 
This suggests that for the reconstruction of 7 the values of the scattering data for small 
^ are very important. This is analogous to observations made in the 2-D case [MS03]. 

5.5. Influence of the support and magnitude 0/7(0;)^/^ — 1 

So far we have used fixed values for d and a, which determine the support and the 
magnitude of 7(0;)^/^ — 1. Figure 5 displays the reconstructions 7'"'^ and of 7(^) 



Direct numerical reconstruction of conductivities in 3-D 



17 




Figure 4. Left: reconstructed Schrodinger potential with truncation of t at i? = 15, 25, 
and 50. Right: reconstructions of 7. 



from t™''(,^) for different choices of support d and magnitude a. Each row corresponds 
to a certain rf-value and each column to a specific «- value. For small support and small 
magnitude we get good reconstructions, but the quality changes dramatically with larger 
amphtude and larger support. Especially 7"''^ does not recover the actual conductivity 
very well for the large amplitude a — 0.9. 

6. Conclusions 

In this work a direct method based on [Nac88] for reconstructing a 3-D conductivity 
distribution from the Dirichlet-to-Neumann map was implemented and tested on noise- 
free data. A linearizing approximation to the scattering transform, denoted t™'" was 
studied and compared to Calderon's reconstruction algorithm. Reconstructions of 
spherically symmetric conductivities in the unit sphere were computed using the 
approximation, Calderon's method, and a scattering transform computed from the 
definition requiring knowledge of the actual Schrodinger potential. The latter case 
served as a benchmark to study the quality of reconstructions for which the actual 
scattering transform is known. It was shown that very accurate reconstructions can 
be obtained from accurate knowledge of the scattering transform. It was found that 
in contrast to the 2-D case, the t""'' approximation is inaccurate near the origin, and 
this results in poor approximations to the magnitude of the conductivity. However, 
the support of 7 — 1 and the boundary value 7 = 1 was well approximated by all three 
methods. Truncating the computed scattering transform in the computations was found 
to have a profound effect on the reconstructed Schrodinger potential q, but the affect 
on the reconstructed conductivity 7 was less dramatic. In summary, it appears that 
the use of the full scattering transform in this method is a promising approach for 3-D 
reconstructions, while linearizations lead to significant inaccuracies in the reconstructed 
amplitudes. 
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Figure 5. Reconstructions of conductivities of varying support and magnitude: each 
row corresponds to a specific support d and each column corresponds to a specific 
magnitude of 7. The dash-dotted curves are the 7""^ reconstructions, the dashed 
curves are the 7*^^^ reconstructions, and soUd curves are the actual conductivities 7. 
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